%%%%%%%%%%%%%%%%%%%%%%Data%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
data = dataset('File','tableC4_data.csv','Delimiter',',');
iter = 1;
warning('off','all');

stations = [data.station_id1, data.station_id2, data.station_id3];

vars = {'num_resp', 'num_vb_parties','turnoutrate_2009','voteshare_vb_2009', 'electors_2014', 'num_cand', 'num_party', 'poll_date1', 'poll_date2', 'poll_date3', 'poll_date4', 'poll_date5', 'state_election', 'lit_pc', 'rural_pc', 'scst_pc'};

%Control mean and SD IPW
fprintf('column 1 \n');
for i = 1:16
    y0 = data.(vars{i});
    idx = find(~isnan(y0) & data.treat==0);
    y = y0(idx);
    num = size(idx);
    wgt = data.wgt_treat(idx);
    treat = data.treat(idx);
    mean_ipw = (sum(wgt.*y))/(sum(wgt));
    sd_ipw = std(y,wgt);
    fprintf('%6.2f\n (%6.2f) \n',mean_ipw, sd_ipw);
end;

fprintf('column 2-4 \n');
for i = 1:16
    y0 = data.(vars{i});
    idx = find(~isnan(y0));
    y = y0(idx);
    num = size(idx);
    wgt = data.wgt_treat(idx);
    treat = data.treat(idx);
    [beta,se,p,r2,F,control_mean,N,C,T] = ipwreg(y,treat,ones(num),wgt,stations(idx,:),2);
    p_ri = ri(iter,beta,se,y,ones(num),stations(idx,:),2,data.prob_treat(idx));
    fprintf('%6.2f (%6.2f) %6.2f %6.2f\n',beta(1),se(1),p(1),p_ri(1));
end;

X = [];
y0 = data.num_resp;
idx = find(~isnan(y0));
num = size(idx);
for i = 1:16
    X = [X data.(vars{i})(idx)];
end;
        
[~,~,~,r2,F0] = ipwreg(data.treat(idx),ones(592,1),X,data.wgt_treat(idx),stations(idx,:),2);
p0 = 1 - fcdf(F0,16,41);

X2 = [ones(592,1) X];
p_ri = f_ri(iter,F0,data.treat(idx),X2,stations(idx,:),data.prob_treat(idx));
fprintf('Joint F-statistic is %6.2f, p-value (standard) is %6.2f, and p(RI) is %6.2f',F0,p0,p_ri);


